clear all
close all


%% chage root to the path of the folder where ReproductionMatlab is located
root="/Users/luigipaciello/Dropbox/RESEARCH/1 PUBLISHED PAPERS/JPE2025/ReproductionMatlab";


%% select experiment
incid_expe = 1;

%set incid_expe=1 to obtain Figures 8-11 & 14, O6 (Appendix); 
%set incid_expe=2 to obtain Figure 12; 
%set incid_expe=3 to obtain Figure 15a & O8 (Apendix),
%set incid_expe=4 to obtain Figure 13; this has to be run after running incid_expe=1
%set incid_expe=5 to obtain Figure 15b,
%all produce Figures 3-6 & O2-O3 (Appendix)



%% Solve the experiment
str1="/calibration_data/";
path_load=append(root,str1);
addpath(path_load)
str2="/results/";
path_save=append(root,str2);
addpath(path_save)
str3="/calibration_Table2/";
path_calibration=append(root,str3);
addpath(path_calibration)
str4="/codes/";
path_codes=append(root,str4);
addpath(path_codes)

%load parameters
calib         = 1;
get_params; 

db   =0.025; %this has to be small enough so to get a good estinate of bbar
dz   =0.1;
z_max=15;
dt   =1;
T_y  =100;

b_vec_mat_l = [-4.6:db:log(1+db),db];
z_vec       = log(max(0.05,dz)):dz:log(z_max);
db_mat      = db;
N_z         = length(z_vec);
N_b         = length(b_vec_mat_l);
Ntot        = N_z*N_b;
get_matrices;

%% FIGURES 3-6 & O2-O3 (Appendix)
%get initial steady state equilibrium in the south
get_firm_policy_loopVA;

%get initial steady state equilibrium in the north
inde_dilution=1;policy_expe=0;
get_equilibrium_north;

%%  Otpimal Subsidy
psi    =-log(.5)/10;
T_sim = T_y/dt;

%matrices for decomposition of different channels 
dilution_vec      =[1 1 1 1 1 0];
l_reallocation_vec=[0 0 1 1 0 1];
demand_ext_vec    =[0 1 0 1 0 1];
shiftdist_vec     =[1 1 1 1 0 1];


W_mat=[];
if incid_expe == 1
    %% Figures 8-11 & 14, O6 (Appendix)
    inde_dilution  = 1;
    l_reallocation = 1;
    demand_ext     = 1;
    shiftdist      = 1;
    k_vec   = 1.0:-0.05:.5;incid_vec=ones(length(k_vec),1);
    tau_vec   = [0.0:0.1:.7,.75,.76,.77,.78,.79,.8,.85];
    ip2_sub=find(k_vec<=0.75,1);ip1_sub=find(tau_vec==0,1);

    get_equilibrium_south;
    W_mat=welf_disc_all;
    k_vec_base=k_vec;
    W_mat_base=W_mat;
    inde_mix = find(tau_vec==0);
    [mm,inde_size_benefit]=max(W_mat(inde_mix,:));
    if min(inde_solved(:,inde_size_benefit))==0,warning('no solution somewhere'),end
    if min(inde_solved(inde_mix,:))==0,warning('no solution somewhere'),end

    get_figures_ss;get_figures_transition;

elseif  incid_expe == 2 
    %% FIGURE 12 & O7 (Appendix)
    k_vec   = 1.0:-0.05:.50;incid_vec=ones(length(k_vec),1);
    tau_vec   = 0:0.1:.8;
    ip2_sub=100;ip1_sub=100;
    W_mat=zeros(length(tau_vec),length(k_vec),6);
    for ie=1:6
        inde_dilution  = dilution_vec(ie);
        l_reallocation = l_reallocation_vec(ie);
        demand_ext     = demand_ext_vec(ie);
        shiftdist      = shiftdist_vec(ie);
        if ie==4 || ie==6
            tau_vec   = 0:0.1:.8;
        else
            tau_vec   = 0;
        end

        get_equilibrium_south;
        if ie==4 || ie==6
            W_mat(:,:,ie)=welf_disc_all;
        else
            W_mat(1,:,ie)=welf_disc_all;
        end
    end
    get_figures_decomposition;

elseif incid_expe == 3
    %% FIGURE O8 (Apendix) & simulations for O9
    incid_ISS      = 0.2;
    inde_dilution  = 1;
    l_reallocation = 1;
    demand_ext     = 1;
    shiftdist      = 1;
    k_vec          = 1.0:-0.05:.5;incid_vec=[1,incid_ISS*ones(1,length(k_vec)-1)];
    tau_vec        = [0,sub_time];
    ip2_sub        = 100;ip1_sub=100;

    get_equilibrium_south;
    W_mat_Iss_size=welf_disc_all;

    inde_sub = find(k_vec<=1-sub_size,1);
    mass_entry_s = [entry_t_s(end,1,1)*ones(1,7),entry_t_s(1:4,2,inde_sub)'];
    mass_entry_n = [entry_t_n(end,1,1)*ones(1,7),entry_t_n(1:4,2,inde_sub)'];

    entry_for_simulations.mass_entry_s = mass_entry_s;
    entry_for_simulations.mass_entry_n = mass_entry_n;
    str="/entry_for_simulations";
    str_save=append(path_save,str);
    save(str_save, 'entry_for_simulations');

    c_ss=cons_ss;
    inde_solved_tau=find(inde_solved(1,:)==1);
    inde_solved_k5=find(inde_solved(2,:)==1);

    figure(100)
    plot(100*(1-k_vec(inde_solved_tau)),100*1/c_ss*(W_mat_Iss_size(1,inde_solved_tau)-W_mat_Iss_size(1,1)),'-b','Linewidth',6),hold on
    plot(100*(1-k_vec(inde_solved_k5)),100*1/c_ss*(W_mat_Iss_size(2,inde_solved_k5)-W_mat_Iss_size(1,1)),'-.r','Linewidth',6),hold off
    xline(100*sub_size,':k','ISS','LineWidth',2,'FontSize',18,'Interpreter','Latex','LabelVerticalAlignment','bottom')
    %yline(0,'--k','Status quo','LineWidth',2,'Interpreter','Latex','Fontsize',18,'LabelHorizontalAlignment','right')
    xlabel('$\vartheta+\tau$, \%','Fontsize',18,'Interpreter','Latex')
    ylabel('Consumption gains, in \% dev. from  status quo','Fontsize',18,'Interpreter','Latex')
    axis([0 30 0 0.25])
    grid on
    le=legend('${\bf{i}} = 0.2$, $\tilde\tau=0.00$','${\bf{i}} = 0.2$, $\tilde\tau=0.79$');
    set(le,'Interpreter','Latex','Fontsize',16,'Location','NorthWest');
    str1="/figures/figure_wefare_incidence_size.pdf";
    str_save1=append(root,str1);
    settt
    print(gcf,'-dpdf',str_save1)
    k_vec_incid= k_vec(inde_solved_tau);
    W_incid     = 1/c_ss*W_mat_Iss_size(1,inde_solved_tau);
  
    %% variation over incidence  figure 15
    incid_vec=[1,1:-0.1:0.1];k_vec = [1,(1-sub_size)*ones(1,length(incid_vec)-1)];
    tau_vec   = [0,sub_time];
    ip2_sub=100;ip1_sub=100;
    debt0t_s=[];debt0_s=[];
    get_equilibrium_south;
    W_mat_Iss_incid=welf_disc_all;

    figure(3)
    plot(100*(incid_vec([inde_solved_k5(3:end)])),debt0t_s(2,[inde_solved_k5(3:end)]),'-b','Linewidth',6),hold on
    plot(100*(incid_vec([inde_solved_k5(3:end)])),debt0_s(2,[inde_solved_k5(3:end)]),'--r','Linewidth',6),hold off
    xline(20,':k','ISS','LineWidth',2,'FontSize',18,'Interpreter','Latex','LabelVerticalAlignment','bottom')
    yline(debt0t_s(1,end),'--k','Status quo','LineWidth',2,'Interpreter','Latex','Fontsize',18,'LabelHorizontalAlignment','right')
    %xlim([10 100 ])
    %ylim([7 18 ])
    grid on
    xlabel('${\bf{i}}$, \%','Fontsize',18,'Interpreter','Latex')
    ylabel('Debt at entry: $b_0$','Fontsize',18,'Interpreter','Latex')
    le=legend('Unsubsidized','Subsidized');
    set(le,'Interpreter','Latex','Fontsize',16,'Location','SouthEast');
    str1="/figures/figure_debt0_tau50.pdf";
    str_save1=append(root,str1);
    settt
    print(gcf,'-dpdf',str_save1)

    %% simulated data for Table 6
    incid_vec=[1,0.2];k_vec = [1,1-sub_size];
    tau_vec   = sub_time;dt=1;
    ip2_sub=100;ip1_sub=100;
    debt0t_s=[];debt0_s=[];
    get_equilibrium_south;if min(min(inde_solved))==0,error('no solution'),end

    % Simulate
    n_years = 11; dt=0.1;
    A_vec_t_s = zeros(n_years/dt,1);
    A_vec_t_n = zeros(n_years/dt,1);
    A_vec_t_s(1:7/dt)=A_s_ben;
    A_vec_t_n(1:7/dt)=A_n_ben;
    A_vec_t_s(7/dt+1:end)=A_s_new;
    A_vec_t_n(7/dt+1:end)=A_n_new;
    w_t_s = zeros(n_years/dt,1);
    w_t_n = zeros(n_years/dt,1);

    time=0;l_s=l_s_ss;l_n=l_n_ss;t_sim=0;
    while time<11-dt/2
        time=time+dt;t_sim=t_sim+1;
        Y_new          = (wei_s*l_s*A_s_new^(1/(1-nu)) + wei_n*l_n*A_n_new^(1/(1-nu)))^((nu-1)/(nu-2));
        w_t_s(t_sim)   = (nu-1)/nu*(A_vec_t_s(t_sim)/Y_new)^(1/(1-nu));
        w_t_n(t_sim)   = (nu-1)/nu*(A_vec_t_n(t_sim)/Y_new)^(1/(1-nu));
        l_s_new  = l_s_ss + dl_s*(1-exp(-psi*time));
        l_n_new  = (1-l_s_new*wei_s)/wei_n;
        l_n      = l_n_new;
        l_s      = l_s_new;
    end

    n_sim_base  =100000;

    %south
    A_iss=A_s_new;A_ss=A_s_ben;
    b_bar_iss=ppval(pp_bar,A_iss);b_bar_ss=ppval(pp_bar,A_ss);
    q_hig=q_hig_s;q_low=q_low_s;q_mid=q_mid_s;
    b_up_unfund=b_up2;b_md_unfund=b_md2;b_dn_unfund=b_dn2;
    b_up_fund=b_up;b_md_fund=b_md;b_dn_fund=b_dn;
    w_t=w_t_s;
    entry_scaled = mass_entry_s/mass_entry_s(1);

    get_sim_regression;
    data_for_regression_s = data_mat;
    data_regression.data_for_regression_s = data_for_regression_s;

    %north
    A_iss=A_n_new;A_ss=A_n;
    b_bar_iss=b_bar_n;b_bar_ss=b_bar_n;
    q_hig=q_hig_n;q_low=q_low_n;q_mid=q_mid_n;
    b_up_unfund=b0_hig_n;b_md_unfund=b0_mid_n;b_dn_unfund=b0_low_n;
    b_up_fund=b0_hig_n;b_md_fund=b0_mid_n;b_dn_fund=b0_low_n;
    b_up_ss=b0_hig_n;b_md_ss=b0_mid_n;b_dn_ss=b0_low_n;
    w_t=w_t_n;
    entry_scaled = mass_entry_n/mass_entry_n(1);

    get_sim_regression;
    data_for_regression_n = data_mat;


    data_regression.data_for_regression_n = data_for_regression_n;
    str="/data_regression";
    str_save=append(path_save,str);
    save(str_save, 'data_regression');


    str="/results/data_n.txt";
    str_save=append(root,str);
    csvwrite(str_save,data_for_regression_n)
    str="/results/data_s.txt";
    str_save=append(root,str);
    csvwrite(str_save,data_for_regression_s)

    str="/results/mass_entry_n.txt";
    str_save=append(root,str);
    csvwrite(str_save,[mass_entry_n(1:11)',(2010:2020)'])
    str="/results/mass_entry_s.txt";
    str_save=append(root,str);
    csvwrite(str_save,[mass_entry_s(1:11)',(2010:2020)'])

elseif incid_expe == 4 
    %% Figure 13
    get_optimal_subsidy_dixit;

elseif incid_expe == 5
    %% Figure 15b
    get_simulation_incidence;
end